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Simulations of single and two-component galaxy 
decompositions for spectroscopically selected galaxies from 
the Sloan Digital Sky Survey 
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ABSTRACT 

We present simulations used to test the two dimensional decompositions of SDSS 
galaxies utilizing the fitting routine GALFIT and analysis pipeline PyMorph. Analy- 
sis showing the bias and scatter of the recovered parameters is presented for multiple 
combinations of simulated galaxy models and fit types. The simulations show that 
accurate measurement of single Sersic models is well constrained when using SDSS- 
quality data. Two component fits present a less robust image decomposition. Galaxies 
observed at higher redshift by Hubble are also simulated. We examine the bias created 
when fitting incorrect models to galaxies. Fitting a two-component Sersic + Exponen- 
tial model to what is really just a single Sersic results in a noisier recovery of the 
input parameters, but these are not biased; fitting a single Sersic to what is truly a 
two-component system results in significant biases. These biases, for total magnitude 
and halflight radius in particular, should be useful in correcting other automatic fitting 
routines. 

Key words: galaxies: structural parameters - galaxies: fundamental parameters - 
galaxies: catalogs - methods: numerical - galaxies: evolution 



• 1 INTRODUCTION 



Analytical descriptions of photometric galaxy structure are 
useful data for evaluating models of galaxy dynamics and 
evolution (e.g., Shankar et al. 2012). Common analytical 
models include the r 1//4 and the r 1//n models originally pre- 
sented by de Vaucouleurs (1948) and Sersic (1963). Model 
development was driven by observations of elliptical galax- 
ies in the case of de Vacouleurs while Sersic observed a con- 
tinuum of early to late-type galaxies. Sersic attempted to 
encapsulate the change in profile shape through use of the 
Sersic index, n, contained in his model. Subsequent empirical 
study suggests that bulges and elliptical galaxies are better 
described by de Vacouleurs profiles or Sersic profiles with 
Sersic index n ~ 4. Disks and late-type spirals are best de- 
scribed by exponential profiles, or Sersic profiles with Sersic 
index n « 1 (Freeman 1970). More recent work has shown 
that the relationship between Sersic index and the photo- 
metric or kinematic components of a galaxy is more compli- 
cated. Following Kent (1985), many analyses simultaneously 
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fit a second component in order to better accommodate the 
qualitative differences of bulges and disks in galaxies. 

There has been recent work on improving photometric 
decomposition of galaxies (Simard et al. 2011; Kelvin et al. 
2012; Lackner & Gunn 2012) with particular interest on the 
applicability to evolutionary models. However, different fit- 
ting processes continue to produce results which differ qual- 
itatively and quantitatively. Systematics in fitting continue 
to be of concern, and the reliability of two-component de- 
compositions in cases of low to moderate signal-to-noise are 
often viewed with some skepticism. In order to quantify the 
systematics and robustness of the fits presented in Meert, 
Vikram & Bernardi (2013), hereafter referred to as M2013, 
we generate simulations of single and two-component galax- 
ies and fit them using the same PyMorph pipeline (Vikram 
et al. 2010) used for the photometric decompositions pre- 
sented in M2013. The M2013 catalog has already been used 
in Bernardi et al. (2012) to study systemmatics in the 
size-luminosity relation, in Shankar et al. (2012) to study 
size-evolution of spheroids, and in Huertas-Company et al. 
(2012) to study the environmental dependence of the mass- 
size relation of early- type galaxies. 

Following several detailed studies which have used sim- 
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ulations to test the robustness of different fitting algorithms 
(e. g., Haussler et al. 2007), the main goal of this paper is to 
use simulations to test the PyMorph pipeline software. We 
present more general tests than Haussler et al. (2007). Our 
simulations are specifically applicable to the SDSS galaxies 
and are useful for evaluating the decompositions presented 
in M2013. 

A description of the simulation process is presented in 
section [5] This includes construction of a catalog of realistic 
galaxy parameters (section 12. generating galaxy surface 
brightness profiles based on these parameters (section |2.2[) ; 
generating sky and noise (sections 12.31 12. 4|) ; and including 
seeing effects in the final image. The completed simulations 
are run through the fitting pipeline and the fits are ana- 
lyzed in section [3] We examine the recovery of structural 
parameters in noise-free images (section I3.1|l and parame- 
ter recovery in realistic observing conditions including both 
neighboring sources and the effects of incorrect PSF estima- 
tion (section 13. 2[) . The effect of changing the fitted cutout 
size (section l3.3[) and the effect of incorrect background esti- 
mation (section l3.4p are also examined. PyMorph is further 
tested by inserting simulated galaxies into real SDSS im- 
ages to test the dependence on local density (|3.5|) . Finally, 
we briefly examine the effects of further lowering the signal- 
to-noise to levels found in the GEMS survey (Haussler et al. 
2007), which used the Hubble Telescope Advanced Camera 
for Surveys (section l3.tif) . In section [3] we discuss the overall 
scatter in our fits and the implications of the simulations. 
Finally, in section [5] we provide concluding remarks. 

We generate single component Sersic galaxy mod- 
els (hereafter referred to as Ser) and two forms of two- 
component galaxy models: one is a linear combination of 
de Vacouleurs and an exponential profile (DevExp) and and 
the other replaces the de Vacouleurs with a Sersic profile 
(SerExp) . A good overview of the Sersic profile used through- 
out this paper is presented in Graham & Driver (2005). 
Throughout the paper, a ACDM cosmology is assumed with 
(h,£l m ,QA) = (0.7,0.28,0.72) when necessary. 

2 CREATING THE SIMULATIONS 

2.1 Selecting the simulation catalog 

We create a set of simulated galaxies using fits from the 
the photometric decompositions presented in M2013. These 
galaxy parameters represent the r-band image decomposi- 
tions of a complete sample of the SDSS spectroscopic cat- 
alog containing all galaxies with spectroscopic information 
in SDSS DR7 (Abazajian et al. 2009). Cuts on the surface 
brightness to ensure completeness and on redshift to pre- 
vent contamination of the estimated distance due to pecu- 
liar velocity are imposed. A brightest magnitude cut is also 
applied to remove galaxies that are too close-by and too well 
resolved to be treated in this simplified manner. The sam- 
ple and the cuts used are similar to the catalog presented 
in Simard et al. (2011); see M2013 for more information on 
the sample of galaxies and the cuts imposed. 

The fit parameters of the photometric decompositions 
were determined using the PyMorph pipeline (Vikram et al. 
2010) which uses SExtractor (Bertin & Arnouts 1996) to 
identify neighboring sources, generates masks, estimates ini- 
tial parameters and fits galaxies using the GALFIT (Peng 



et al. 2002) 2D fitting routine. A number of other non- 
parametric measurements are also reported. 

From the spectroscopic sample, three subsamples of 
galaxies that are representative of the SDSS spectroscopic 
catalog and produce physically meaningful photometric de- 
compositions were selected. In order to ensure that the 
galaxies are representative of the full catalog, we exam- 
ined the distributions of basic observational parameters of 
SDSS galaxies (surface brightness, redshift, apparent Pet- 
rosian magnitude, Petrosian half-light radius, and absolute 
magnitude) . 

Some restrictions on fit parameters are necessary to en- 
sure that outliers are removed from the parameter space 
used to generate the simulations. Galaxies that do not sat- 
isfy these basic cuts are removed to ensure that the param- 
eters used to generate the images are not unphysical. The 
cuts do not significantly bias our galaxy distribution as is 
shown in Figure [1] The cuts include: 

(i) Any Sersic components must have Sersic index less 
than 8. 

(ii) Half-light radius of any Sersic component must be less 
than 40 kpc. 

(iii) In the two-component fits, the ratio of the bulge 
halflight radius to disk scale radius should be less than 1, or 
the galaxy should be bulge dominated (B/T > 0.5). 

Conditions (i) and (ii) were chosen to prevent selection of 
Ser models with extended profiles that are likely the result 
of incorrect sky estimation during the fitting process. Condi- 
tion (iii) is used to ensure that any disk dominated galaxies 
have a bulge component that is smaller than the disk. 

After enforcing the cuts on the sample of fitted galax- 
ies, 10,000 galaxies for each of the Ser, DevExp, and SerExp 
samples were selected at random. The fitted parameters of 
these sample galaxies were used to generate simulated model 
galaxies used in testing the pipeline. 

Figure [T] shows the distributions of surface brightness, 
redshift, extinction-corrected r-band Petrosian magnitude, 
r-band Petrosian halflight radius, and absolute magnitude of 
all SDSS spectroscopic galaxies (in black) and our simulation 
samples: Ser (red), DevExp (green), and SerExp (blue). The 
distribution of simulated galaxies reproduces the observed 
distribution for all three samples for each observational pa- 
rameter as verified by a KS 2-sample test. 

Figure [T] also presents the signal-to-noise (S/N) of the 
simulation samples as compared to the parent distribution. 
The S /N of the images is a limiting factor in the fitting pro- 
cess, so care must be taken to ensure that the S/N is not 
artificially increased in the simulations when compared to 
true SDSS galaxies. This S/N is calculated using the r-band 
Petrosian magnitude and r-band Petrosian halflight radius. 
Petrosian quantities are used to make a fairer comparison 
among all the samples. Because the Petrosian quantities are 
non-parametric, they avoid the complications that arise in 
assessing the possible biases introduced during fitting. Any 
differences in S/N are not large enough to significantly bias 
the distributions as verified by the KS 2-sample test. There- 
fore, we conclude that our samples are fair representations of 
the underlying distribution of SDSS spectroscopic galaxies. 

An unbiased selection in the previously mentioned pa- 
rameters is not sufficient to guarantee fair sampling of the 
S/N with respect to magnitude, nor does it prevent ficti- 
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band Petrosian halflight radius, [(e)| V m ax-weightcd luminosity function, and | (f ) | signal- to- noise distribution of the simulated sample 
of galaxies drawn from the SDSS spectroscopic galaxy sample of DR7. The distribution of all SDSS spectroscopic galaxies is shown in 
black. Distributions of the Ser, DevExp, and SerExp simulations are shown in red, green, and blue, respectively. Bin counts are normalized 
to integrate to 1. The distributions of the simulated galaxies are representative of the full sample fitted in M2013 and are appropriate 
to compare to the SDSS spectroscopic sample as verified by a Kolmogorov-Smirnov 2-sample test. The signal-to-noise (S/N) will be 
discussed further in section \1. 41 In calculating this S/N, we use the measurement of sky provided by the PyMorph pipeline rather than 
SDSS to identify and separate target counts from sky counts. PyMorph sky estimation is shown to be more accurate than the SDSS 
estimation provided in the DR7 catalog. 



tious correlations among multiple fit parameters. In fact, 
correlations among fit parameters are to be expected if the 
PyMorph pipeline is robustly measuring properties of the 
target galaxies (there are many obvious correlations that ex- 
ist among parameters, as evidenced by the numerous scaling 
relations that currently exist). However, it is difficult, and 
largely unnecessary, to examine every possible relationship 
thereby ensuring that all correlations are true reflections of 
the underlying data. 

Two important relationships are the halflight radius as 
a function of apparent magnitude and the S/N as a function 
of the apparent magnitude. The two relationships are neces- 
sary to ensure the appropriateness of the simulation because 
systematic differences in radius will lead to systematic vari- 
ation in the S/N of the sample. Scatter in the recovered 
fitting parameters as a function of magnitude is particularly 
useful, since it is relatively easy to estimate total magnitude 
on image data. Therefore, the S/N as a function of apparent 
magnitude should appropriately reflect that of the parent 
sample from SDSS. 

Figure [5] presents the halflight radius versus apparent 
magnitude, and Figure [3] presents the S/N versus apparent 
magnitude. The simulated samples are shown in red, green, 
and blue corresponding to the Ser, DevExp, and SerExp mod- 
els. The underlying SDSS parent distribution is shown in 
black. Indeed, Figure [2] shows that the Ser and SerExp mod- 
els are in close agreement with the full SDSS sample. The 
SerExp sample has smaller radii (about 10% smaller) in the 
brightest magnitude bin (14.0 < mtot.r < 14.5). The Ser 
model radii agree across the magnitude range. However, the 
S/N agrees with the full SDSS sample, or is slightly be- 
low that of SDSS. The lower signal-to-noise, although not 
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Figure 2. The distribution of galaxy radii as a function of ap- 
parent magnitude for the parent SDSS sample in black, the Ser 
model in red, the DevExp model in green, and the SerExp model 
in blue. The median in each bin is shown with corresponding 
contours showing the extent of the middle 68% of data (dashed 
lines) and middle 95% of the data (dot-dashed lines). At brighter 
magnitudes, the DevExp model begins to diverge from the orig- 
inal sample toward smaller radi. The SerExp model is in closer 
agreement across the entire magnitude range. 



exactly that of SDSS will not bias the tests toward better 
results, so we deem these samples acceptable for testing. 
The DevExp model tends to smaller radii and higher S/N at 
brighter magnitude. The DevExp model was also analyzed, 
however the results are not discussed in this paper. They 
can be found in Meert et al. (2013). 

Testing the accuracy of the PyMorph fitting routine 
does not necessarily require an unbiased parameter distri- 
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Figure 3. The distribution of galaxy S/N as a function of appar- 
ent magnitude is presented in the same format as Figure [2] At 
brighter magnitudes, the DevExp model begins to diverge from the 
original sample toward higher S/N. The SerExp and Ser models 
are in closer agreement across the entire magnitude range. 



bution. In reality, all that is required is a sample with suf- 
ficient coverage of the parameter space represented by the 
data. The simulations use smooth profiles, simplifications of 
the true galaxies that are observed in SDSS. Examination of 
the results of fitting these simplified models and comparison 
to fits of true observed galaxies can potentially yield useful 
information regarding galaxy structure. In Bernardi et al. 
(2012), the simulations are used together with the decom- 
positions of the SDSS spectroscopic sample to characterize 
the scatter in the size-luminosity relation as well as examine 
possible biases. In order to make these comparisons, an un- 
biased sample is required. The distributions shown in Figure 
[T]show that the simulations are appropriate to use for this 
purpose. 



2.2 Generating the images 

We generate the two-dimensional normalized photon distri- 
butions from the one-dimensional Sersic profiles and the one- 
dimensional exponential profiles of each bulge and disk com- 
ponent. Disk components are only simulated where required, 
as is the case for two-component fits. When multiple com- 
ponents are to be simulated, each component's normalized 
photon distribution is generated separately and combined 
prior to generating the simulated exposure. 

Two-dimensional galaxy profiles are treated as az- 
imuthally symmetric one-dimensional galaxy light profiles 
that are deformed according to an observed ellipticity. The 
galaxy profiles are generated using the structural parameters 
generated from photometric decompositions as described in 
the previous section. Single component galaxy profiles and 
the bulges of two component galaxies are generated accord- 
ing to the Sersic profile 



J(r) = I e exp 



R e 



(1) 



b n = 1.9992n - 0.3271 



where Sersic index (n), half-light radius (-Re), and surface 
brightness at R e (I e ) are selected simultaneously from the 
catalog described in the previous section. 

For the DevExp and SerExp cases, an exponential disk 



(the formula above, with n = 1) is added to the Sersic com- 
ponent to model the disk component of the galaxies. This 
model requires input parameters scale radius (Rd) and cen- 
tral surface brightness (Id). The disk is modeled using 



(2) 



Iex P {t) = hexp — 

ttd 



The disk profile can be recast as equation JT]) by noting that 
the exponential profile is a special case of the general Sersic 
profile. The equivalence is true when 

I d = 7 e exp&„, R e = (b n ) n Rd, and n = 1. (3) 

The one-dimensional profiles are made into two dimen- 
sional profiles by assuming a uniform distribution in angle 
9. 

At this point, there is no explicit dependence in the pro- 
files on the position angle, <j>, of the components. Therefore, 
we begin by treating each component as having the major 
axis coincident with the x-axis. For example, a photon with 
8 — will fall on the semi-major axis of the galaxy while a 
photon with 8 — tt/2 will fall on the semi-minor axis. 

The continuous image in (r, 9) space is then transformed 
into (x, y) space in the plane of the sky using the axis ratio 
(b/a) of the respective component according to the equations 



x = (b/a) 1 ^ 2 r cos 9 



V = (b/a) 



-1/2 



r sin 9 



(4) 



At this point, the component major axis is still coinci- 
dent with the positive x axis. To rotate by the true position 
angle of the component (cpd for the disk or cpi for the bulge), 
we use the standard rotation matrix 



R : 



COS (f> 

sin (f> 



— sin ( 

cos tb 



(5) 



where <j> — 4>t or 4>d depending on which component is cur- 
rently being simulated (bulge or disk component, respec- 
tively) . 

For true bulges, the axis ratio defines an intrinsic ellip- 
ticity of the bulge. More exactly, this ellipticity is believed 
to be a two-dimensional projection of an intrinsically ellip- 
tical object. On the other hand, we treat the exponential 
portions of galaxy profiles as an infinitesimally thin disk. 
Rotationally supported thin disk components are thought 
to be intrinsically round, and so their ellipticity is not inter- 
preted as a true ellipticity but rather a projection effect due 
to the line of sight used to view the disk. For this reason, 
the disk ellipticity is often expressed as an inclination, ipd, 
where 



ipd = cos 1 {b/a) 



(6) 



so that b/a = 1 implies ipd = (the observer is viewing the 
disk face-on) and b/a = implies ipd = 7r/2 (the observer is 
viewing the disk edge-on). 

After rotation, the continuous probability distribution 
represented by the two-dimensional light profile can be in- 
tegrated over pixel areas to define a pixellized probability 
distribution. The details of this integration are mainly unim- 
portant. However, the simulation must take careful account 
of the integration in the central pixels, where the profile 
can vary greatly over a single pixel. Various oversampling 
methods have been devised to properly correct this com- 
mon problem (e.g., Peng et al. 2002; Haussler et al. 2007). 
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Figure 4. The distribution of sky values for data in the 
SDSS CASJOBS catalog. These data are drawn from the sky_r 
CASJOBS parameter and are converted into counts (DN) per 
pixel per standard SDSS image exposure of 54 seconds. We use 
this distribution to determine the sky value used in our sim- 
ulations. As an approximation, we use the mean value of 130 
counts / pixel / exposure. 

The simulations in this paper have been tested to ensure 
that the pixel-by-pixel integration is accurate to w 3% of 
the corresponding Poisson noise in a given pixel. Therefore, 
we treat the simulations as exact calculations of the galaxy 
photon distributions. 

After generating image, the galaxy is numerically con- 
volved with a PSF extracted from SDSS DR7 data using 
read_PSF program distributed by 

2.3 Creating the background 

Two hundred background images, each equal in size to an 
SDSS fpC image, are also simulated for testing purposes. 
These images contain constant background and a randomly 
selected field of galaxies taken from an SDSS fpC image. The 
SDSS catalog provides rudimentary photometric decompo- 
sitions of each star and galaxy. Galaxies are fit with an ex- 
ponential disk and a de Vacouleurs (n = 4) bulge indepen- 
dently. The best fit is reported as a linear combination of the 
two fits using the f racdev parameter to express the ratio of 
the de Vacouleurs model to the total light in the galaxy. 

For the simulated background used in this paper, each 
galaxy is generated using the combined profile of the two fits. 
The de Vacouleurs bulge and exponential disk component 
are separately simulated according to the magnitude, radius, 
ellipticity, and position angle reported in SDSS. Each com- 
ponent is simulated using the method described in section 
12.21 The background galaxy is constructed by adding the 
two components using the f racdev parameter. The galaxy 
is then inserted into the fpC image. Any foreground stars are 
also simulated as point sources and inserted into the image. 

For the background sky counts in the image, we use the 
mean sky of all SDSS observations as given in the SDSS 
photoobj table by the sky_r parameter. The distribution 
of the sky flux is plotted in Figure [4] in units of counts (or 
DN) per pixel per exposure. The median and mean values 

1 read_PSF is part of the readAtlasImages-v5_4_ll package avail- 
able at http : //www. sdss . org/dr7/products/images/read_psf .html 



for a 54 second SDSS exposure are ~125 and «130 counts 
per pixel, respectively. We use the mean value of 130 counts 
per pixel as the background in our simulations. This sky 
background is applied to the entire chip as a constant back- 
ground; no gradient is simulated across the image. Back- 
ground gradients should be approximately constant across 
a single galaxy. This assumption is verified by inserting the 
simulated galaxies into real SDSS fpC images near known 
clusters, where the sky contribution should be high and gra- 
dients more likely. In section [331 we show that there is little 
change in the behavior of the fits in these types of environ- 
ments. 

Previous work has improved the measurements of sky 
background (see Blanton et al. 2011). However, these correc- 
tions are largely focused on areas of large, bright galaxies or 
on making the sky subtraction stable for purposes of tiling 
fpC images together. Since we are only focused on main- 
taining the proper S/N for our simulations, the sky levels 
provided in the SDSS database are sufficient, provided that 
they maintain the correct S/N. We discuss the S/N distribu- 
tion of our simulations and the original SDSS galaxy sample 
in section [2741 below. 

Diffraction spikes and other image artifacts are not 
simulated. However, the SDSS photo pipeline often mis- 
identifies additional phantom sources along any observed 
diffraction spike. These sources may be modeled in our back- 
ground, and so these effects may, in fact, be approximately 
modeled, albeit unintentionally. It is reasonable to expect 
that the diffraction effects should not have a large effect on 
the fitting process, as their elongated straight structure does 
not mimic galaxy structure. The dominant effect produced 
by the bright stars in the field is bias in the background 
estimation in the nearby neighborhood of a star. 

After simulation of the background images, and prior 
to adding noise, each background image is convolved with a 
random SDSS PSF selected from original fpC image upon 
which the individual image is based. 

2.4 Noise 

After generating a target galaxy and inserting it into a back- 
ground, Poisson noise is added using the average inverse gain 
of an SDSS image (4.7 e~/DN) and the average contribution 
of the dark current and read noise, referred to as the "dark 
variance," (1.17 DN 2 ), to determine the standard deviation 
for each pixel. Specifically, 

Fi,j = + bkrd»,j (7) 

is the total flux in pixel (i.e., , the sum of the source 

and background fluxes in the pixel), and 

[f~. 

ffi * = 4 / — — + dark variance (8) 
V gain 

so 



for a single pixel. 

Since the fitting pipeline is dependent on the S/N, it 
is essential that the simulated S/N is comparable to SDSS. 
The distribution of the average S/N per pixel within the 
halflight radius for the simulations and the original galaxies 
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is plotted in Figure 1(f) The S/N distribution of simulations 



and the SDSS spectroscopic galaxies agree as verified by a 
KS 2-sample test, therefore the simulations appropriately 
approximate the S/N of SDSS galaxies contained in M2013. 

2.5 Final processing for fitting 

After constructing the images, we generate an image of the 
(Tij values according to equation [8] This image is supplied 
along with the input image to the pipeline in order to cal- 
culate the x 2 value for the fit. 

Figure [5] shows some examples of simulated galaxies 
throughout the simulation process. This includes the noise- 
less simulated galaxy, the noiseless simulated background, 
the composite image of galaxy and background, and the 
composite image after adding Poisson noise with ffij de- 
fined in equation [8] The final image size used for fitting is 
20 times the Petrosian r-band halflight radius. A discussion 
of this choice of stamp size is presented in section 13.31 



3 TESTING PYMORPH IMAGE 
DECOMPOSITIONS 

PyMorph was used to fit the simulated galaxies described 
above in section [2] In order to describe and quantify the 
amount of bias and scatter in the PyMorph image decom- 
position pipeline when applied to SDSS simulations, several 
tests of increasing complexity were performed. The primary 
focus is on the systematics of the pipeline. A better under- 
standing of these systematic effects is useful in interpreting 
the data presented in M2013. 



3.1 Noiseless images 

As an initial test, the pipeline is applied to galaxies with no 
noise and no neighboring sources or background counts. This 
produces the minimum scatter in the data, serves to verify 
that our simulations are correct, and shows that PyMorph 
is properly functioning. 

The total apparent magnitude, halflight radius, and ad- 
ditional fit parameters recovered by fitting the noiseless im- 
ages of the Ser and SerExp models are presented in Figures 

The plots show the dif- 



6(a) 7(a) 8(a) 9(a) and 10(a) 



ference in simulated and fitted values (fitted value - input 
value) versus the input value. The greyscale shows the den- 
sity of points in each plane with red points showing the me- 
dian value. Blue dashed lines show the regions which contain 
68% of the objects. 



Figures 6(a) and 9(a) show the corresponding fit is well 



constrained (Ser fit with Ser, and SerExp with SerExp) al- 
though the scatter increases somewhat for the components 
of the SerExp fit (see Figure 10(a) |. This is most true as the 



component becomes dim (bulge/disk magnitude approaches 
18) where S/N for the individual component drops more 
rapidly. The magnitude and halflight radius are also well 
constrained when a Ser galaxy is fit with a SerExp profile 
(Figure 7(a) l . However, a SerExp galaxy fit with a Ser profile 



produces large biases in the magnitude and halflight radius 



(Figure 8(a) I 



The total magnitude and halflight radius are well con- 
strained (crtotai mag ~ 0.01 mag and <7 r adius » 0.1 arcsec) 



in cases where the correct model is applied to the simu- 
lated galaxy (i. e., Ser galaxy fit with a Ser model). This is 
not necessarily the case when the wrong model is applied 
(i.e., SerExp galaxy fit with a Ser model). It is obvious 
that a single-component galaxy cannot properly model a 
two-component galaxy in general, and therefore, we would 
expect significant problems in attempting to fit a single com- 
ponent profile to a two-component galaxy. Nevertheless, this 
type of fit is often performed on real data when it is believed 
that the resolution and S/N are such that it is unlikely to 
recover a robust two-component fit. An important observa- 
tion is that the SerExp fit provides the most stable estimate 
of the halflight radius and total magnitude regardless of the 
true simulated galaxy model (Ser, DevExp, or SerExp). The 
additional freedom in the SerExp model and the fact that 
the Ser and DevExp models are special cases of the SerExp 
model would lead us to expect this result. 

One noticeable systematic in the pipeline, one that has 
proven difficult to eliminate, is the underestimate of Sersic 
index at larger Sersic indexes. At Sersic indexes of n m 4, 
we underestimate the Sersic index by less than 1%. However, 
this underestimate increases in the later tests. 

Since no Poisson noise is added to these images, the 
effective S/N is much larger than the S/N of an SDSS im- 
age that the PyMorph pipeline is expected to analyze. The 
scatter apparent in these fits is then a combination noise 
associated with the use of the PyMorph routine (as well as 
the GALFIT routine used by PyMorph) and any parameter 
degeneracies inherent in the models. 

GALFIT uses the Levenberg-Marquardt minimization 
method (Press et al. 1992) to find the minimum of the \ 2 
distribution of the fit. The Levenberg-Marquardt method is 
not a global search algorithm but rather follows the steep- 
est decent to a local minimum. As the parameter space be- 
comes more complicated, GALFIT has more trouble accu- 
rately recovering parameters. Adding components to the fit 
(i.e., going from a one-component to two-component fit or 
going from a fixed Sersic index component to one with a free 
Sersic index) will not only complicate the x 2 surface, making 
convergence less likely, but may introduce true degeneracies 
in the parameter space. 

For instance, the SerExp fit of a galaxy of very late 
type often suffers from over-fitting. The bulge component 
will tend to fit the disk of the galaxy as a second disk com- 
ponent with n » 1. This is obviously an unintended solu- 
tion to the fitting, but one that is equally valid from a x 2 
perspective. In practice, it is difficult to prevent this type 
of convergence without artificially constraining the fitting 
routine. Such constraints are generally discouraged and can 
lead to other negative effects including convergence to a non- 
optimal solution. The best solution to the parameter degen- 
eracy is close examination of any two-component fits in cases 
where n « 1, or B/T RiOor 1-B/Tw 0. 



3.2 The effects of background, neighbor sources, 
and incorrect PSF extraction 

When analyzing real data, it is not possible to extract the 
PSF of the target galaxy to arbitrary accuracy. Interpolation 
is required and generally performed on a network of the 
nearest stars to the target galaxy. We test this effect through 
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background 



DevExp model 



background 




model + background+ noise 



SerExp model 



background 



model + background 
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Figure 5. Examples of simulated galaxies and background shown before and after adding Poisson noise. Top, middle and bottom rows 
show randomly selected sample Ser, DevExp, and SerExp profiles, respectively. From left to right, the columns show the simulated galaxy, 
simulated background, background+galaxy, and final image with noise. 



extraction of a neighboring PSF to be used during fitting in 
place of the PSF used to generate the image. 

The neighbor PSF used in fitting is randomly selected 
from a location within a 200 pixel box surrounding the 
source. This provides approximately even sampling of dis- 
tances from nearly to about 170 pixels in separation from 
the source which corresponds to a separation of ~ to 
« 67.32 arcseconds between the target galaxy and the lo- 
cation used for PSF extraction. This inserts some PSF error 
into the process of fitting as would be expected in the case 
of real data. However, it also retains the similarity between 
the PSF used for simulation and the PSF used for fitting. A 
strong similarity between the two would be expected since 
the PSF generally will not vary greatly over the area of a 
single fpC image. 

The PSF effect is applied to target galaxies randomly 
inserted into the simulated fpC images described in section 
12.31 The simulated fpC images contain sky as well as neigh- 
boring sources. These fits are the closest simulation to ac- 
tual observing conditions that we have analyzed. Therefore 
the fits and the resulting measures of scatter and bias are 
adopted as our fiducial estimates of scatter and bias when 
using the pipeline. 



Figures 6(b) 7(b) 8(b) 9(b) and 10(b) show that 



we recover the input values with marginal scatter. The to- 
tal magnitude and halflight radius remain well constrained 
(ototaimag ~ 0.05 mag and <7 ra di us ~ 0.25 arcsec) in cases 
where the correct model is fit to the simulated galaxy. How- 
ever, this scatter becomes much larger when the wrong 
model is fit. The underestimate of the Sersic index, par- 
ticularly at large values, persists. 

Further examination of the two component fits show 
that the pipeline has difficulty extracting dim components 
(bulge or disk magnitude dimmer than w 18). In these 
ranges, the components are observed at lower S/N and the 
pipeline looses sensitivity to the model parameters. The 
SerExp fit shows an underestimate of Sersic index, which 
is even stronger than in the single component case, and an 
underestimate of bulge radius. However, the disk parameters 
remain unbiased with an increase in scatter of the model pa- 
rameters. Most likely, this stability is the result of the fixed 
Sersic index of the bulge (n = 1) which fixes the shape of 
the disk light profile. 

In general, the SerExp fits are problematic and require 
much care when analyzing individual components. However, 
as we have already shown, total magnitude and halflight 
radius are still tightly constrained. 

Table [T] summarizes the bias and scatter in the fits; they 
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Figure 6. The simulated and recovered apparent magnitude, halfiight radius, and Sersic index for a Ser galaxy fit with a Ser model in 
four cases: [(a)] the image prior to adding Poisson noise, [(b)] our fiducial case containing simulated sky and neighboring sources, |(c)| the 
simulated galaxies inserted into real SDSS fpC images, and |(d)| thc simulated galaxies inserted into real SDSS fpC images near clusters. 
Over-plotted are the bias (red points) in the fitted values. All plots show the 68% (dashed line) scatter in blue. The density of points is 
plotted in gray-scale. The Sersic index shows increasing underestimate up to si0.5 (or «8%) at the largest Sersic indexes. 



exhibit trends with both the input value of the parameter 
and the input magnitude of the galaxy. This behavior is not 
properly encapsulated in the overall measure of bias, so these 
values are useful only as an example of the relative scale of 
bias and scatter for each parameter. 



Errors can be correlated across many fit parameters, so 
we also calculate a correlation matrix for the parameter er- 
rors. Figure [TT] shows an example of the correlation matrix 
for the simulated SerExp galaxies fit with a SerExp model. 
We see the expected strong correlations between bulge-to- 
light ratio and the bulge and disk magnitudes as well as the 
correlation among the radii of the bulge component with 
the Sersic index. While the correlation matrix suggests that 
there is little correlation between sky estimation error and 
the fitted parameters, we will show later that there is indeed 
a strong correlation in model errors with sky estimation er- 
ror. 



3.3 Effect of cutout size 

We select postage stamp cutouts for use in fitting. It is im- 
portant to select a cutout size that does not significantly 
bias the fits produced by PyMorph. The most important 
consideration is to provide enough sky pixels to allow the 
PyMorph program to properly determine the sky level in the 
images. Reducing cutout size may cause overestimation of 
background and corresponding errors in the other fit param- 
eters. However, we use the PyMorph pipeline and GALFIT 
to fit a constant background to the galaxy image. Since a 
larger image could make sky gradients more significant, this 
could bias the fits when a larger cutout is used. We seek to 
minimize error when estimating the sky level without intro- 
ducing a gradient term and further complicating the fitting 
process. 

To test for optimal cutout size, we fit simulated galax- 
ies with cutout sizes between 10 and 25 Petrosian half-light 
radii. We plot the average difference between simulated and 
measured fit parameters below. In Figure [12] we present the 
error and 1-a scatter in the error on the total magnitude, 
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Figure 7. The simulated and recovered apparent magnitude and halflight radius for a Ser galaxy fit with a SerExp model in four cases: 
| (a) | the image prior to adding Poisson noise, [(b)] ou r fiducial case containing simulated sky and neighboring sources, |(c)| the simulated 
galaxies inserted into real SDSS fpC images, and |(d)| the simulated galaxies inserted into real SDSS fpC images near clusters. Over-plotted 
are the bias (red points) in the fitted values. All plots show the 68% (dashed line) scatter in blue. The density of points is plotted in 
gray-scale. 
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Figure 8. The simulated and recovered apparent magnitude and halflight radius for a SerExp galaxy fit with a Ser model in four cases: 
| (a) | the image prior to adding Poisson noise, [(b)] our fiducial case containing simulated sky and neighboring sources, |(c)| the simulated 
galaxies inserted into real SDSS fpC images, and |(d)] the simulated galaxies inserted into real SDSS fpC images near clusters. Over- 
plotted are the bias (red points) in the fitted values. All plots show the 68% (dashed line) scatter in blue. The density of points is plotted 
in gray-scale. The inability of the Ser profile to accurately model a SerExp galaxy is clearly evident. Errors in magnitude and halflight 
radius are correlated and the error in radius is largely driven by errors in the largest, brightest objects. However, systematic errors occur 
even at the dimmer magnitudes. Ser fits tend toward recovering larger, brighter objects when applied to a true two component galaxy. 
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Figure 9. The simulated and recovered apparent magnitude, halflight radius, and bulge-to-total light ratio for a SerExp galaxy fit with a 
SerExp model in four cases: [(a)] the image prior to adding Poisson noise, [(b)] our fiducial case containing simulated sky and neighboring 
sources, |(c)| the simulated galaxies inserted into real SDSS fpC images, and |(d)] the simulated galaxies inserted into real SDSS fpC 
images near clusters. Over-plotted are the bias (red points) in the fitted values. All plots show the 68% (dashed line) scatter in blue. 
The density of points is plotted in gray-scale. The apparent trend in B/T from overcstimation at lower B/T values to underestimation 
at higher B/T values is largely due to the boundaries on the parameter space forcing the scatter to be asymmetric (e.g., a galaxy with 
true B/T= f cannot be estimated to have B/T> 1). 
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fitted model 








Ser 




SerExp 






bias 


1-CT 


bias 


1-cr 


Ser 


m to t [mag] 
t m [arcsec] 

sky [%] 
Sersic Index 


0.00 ±0.09 
-0.01 ±0.68 
-0.05 ±0.14 [ 
-0.08 ±0.57 


-0.02,0.02] 
-0.14,0.01] 
-0.11,-0.05] 
-0.44,0.01] 


-0.02 ±0.18 
0.11 ± 1.53 
-0.05 ±0.14 


[-0.11,0.02] 
[-0.11,0.04] 
[-0.11, -0.05] 


SerExp 


mtot [mag] 
r hi [arcsec] 
sky [%] 
B/T 
m 6lli9e [mag] 
m diak [mag] 
tbulge [arcsec] 
r<Msfc [arcsec] 
Sersic Index 


-0.07 ±0.18 [ 
0.49 ± 1.86 
-0.08 ±0.15 [ 


-0.22,-0.01] 
-0.02,0.30] 
-0.16,-0.07] 


-0.02 ±0.15 
0.07 ± 1.17 

-0.06 ±0.13 
0.00 ±0.15 

-0.14 ±0.71 

-0.04 ±0.50 
0.08 ±0.97 
0.07 ±0.82 
0.06 ± 1.98 


[-0.09,0.01] 
[-0.07,0.04] 
[-0.11,-0.05] 
[-0.07,0.03] 
[-0.73,0.06] 
[-0.41,0.05] 
[-0.27,0.11] 
[-0.14,0.08] 
[-0.90,0.14] 



Table 1. The bias and scatter of the fitted parameters of the simulated images with background and PSF effects. These values are 
provided for illustrative purposes only. There is much underlying structure in the errors when compared to their respective input values 
or the magnitude of the component. 
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Figure 10. The simulated and recovered component parameters for a SerExp galaxy fit with a SerExp model in four cases: [(a)] thc image 
prior to adding Poisson noise, [(b)] ou r fiducial case containing simulated sky and neighboring sources, |(c)| the simulated galaxies inserted 
into real SDSS fpC images, and |(d)] the simulated galaxies inserted into real SDSS fpC images near clusters. Over-plotted are the bias 
(red points) in the fitted values. All plots show the 68% (dashed line) scatter in blue. The density of points is plotted in gray-scale. 
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Figure 11. The correlation matrix for a simulated SerExp galaxy 
fit with a SerExp model. 

halflight radius and sky (showing SExtractor sky in blue and 
our estimates in red) as a function of cutout size. Smaller 
sizes clearly bias sky estimates made by SExtractor, but 
only minor improvement in the scatter of any parameters 
is achieved by using cutout sizes above 20 halflight radii. 
Since we use SExtractor sky as a starting point for our fit- 
ting, we choose a size of 20 halflight radii for our images. The 
sky estimates of SExtractor improve substantially. However, 
GALFIT sky estimation is stable over these sizes. Because 
GALFIT sky estimation is largely independent of the initial 
starting SExtractor value (which we would expect if we are 
truly finding the best fit to the galaxy), it is likely the case 
that cutout sizes smaller than even 10 halflight radii could 
be used for analysis. 

Additional plots of other parameters are omitted in this 
section. The other fitted parameters show little or no sensi- 
tivity to cutout size in the range of cutout sizes used. How- 
ever, as previously discussed, the bias and scatter may not 
be equally affected across all model parameters. The effects 
may be concentrated in a small part of the parameter space. 

3.4 The effect of incorrect sky estimation 

Estimation of the sky in the vicinity of the target galaxy has 
a high level of uncertainty. Indeed, accurate sky determina- 
tion is likely not even a solvable problem as discussed briefly 
in Blanton et al. (2011). To determine the bias introduced 
by our sky estimation, we have tested our fitting pipeline 
in cases of both underestimation and overestimation of the 
sky. We fix the sky at the simulated sky level, as well as at 
simulated sky level ±0.5% and ±1.0%. These ranges were 
chosen to represent the range of differences between our sky 
estimations and those provided in the CASJOBS database 
for the SDSS spectroscopic sample. 

Figure[l3]shows a comparison of sky estimates using Py- 
Morph to those provided from the SDSS photometric data 



Figure 13. The percent difference between the sky estimate of 
PyMorph for SDSS galaxies and the sky estimated by the SDSS 
photometric pipeline for those same galaxies. The normalized dis- 
tribution of differences is shown in bins of 0.1%. A negative differ- 
ence indicates that the sky measured by PyMorph is lower than 
that reported by SDSS. The vertical red solid line indicates the 
median of the distribution. The red dashed, dot-dashed, and dot- 
ted lines indicate the 68-95-99% ranges of the data, respectively. 
The 95% range of sky values is approximately between ±1% dif- 
ference, so we adopt this as the range used to test the effects of 
imporper sky estimation. 



pipeline. This comparison is performed on data from the 
catalog presented in M2013. The Figure shows the normal- 
ized distribution of differences in sky estimation in bins of 
0.1%. A negative difference indicates that the sky measured 
by PyMorph is lower than that reported by SDSS. The ver- 
tical red solid line indicates the median of the distribution. 
The red dashed, dot-dashed, and dotted lines indicate the 
68-95-99% ranges of the data, respectively. The 95% range 
of sky values is approximately between ±1% difference. For 
the test, we adopt this range as the range to test for sky 
variation. 

The results of incorrectly estimating the sky are shown 
in Figures [14] and 1151 In red, we show the results of fit- 
ting galaxies using the standard PyMorph pipeline, treating 
sky level as a free parameter in the fit. PyMorph system- 
atically underestimates the sky at the 0.1% level. However, 
the scatter is very tight as indicated by the vertical dashed 
red lines. In black we have plotted the fitting results at fixed 
sky levels of the correct value and ±0.5% and ±1.0%. Er- 
rors approaching 0.5% clearly introduce a large bias in the 
fits. The 0.5% level is an important level because it is the 
approximate level of overestimation shown in the preceding 
section (section |3.3[) found by SExtractor. 

Note the asymmetry of the effects of incorrect sky es- 
timation on fitting parameters. In particular, an underesti- 
mate of sky is much more detrimental to the fit than the cor- 
responding overestimate, because it can have a substantial 
effect on fitted halflight radii and Sersic index. Our pipeline 
tends to underestimate the sky at the «0.1% level which 
would potentially bias our fits. However, this bias is much 
smaller than the bias associated with using the SExtrac- 
tor sky estimate. SerExp disk components are remarkably 
robust to the errors, while bulge paramteres suffer greatly, 
especially when the sky is underestimated. 

Additionally, Figure [13] shows that PyMorph consis- 
tently estimates the sky at a level lower than that of the 
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Figure 12. The mean difference of the total magnitude (left column), PSF-corrected halflight radius (center column), and sky estimation 
(right column) as a function of cutout size for simulated SerExp galaxies fitted with a SerExp model. Other simulated models behave 
similarly. For sky estimation, the sky measured by GALFIT is plotted in red. SExtractor sky measurements are plotted in blue for 
reference. One-cr scatter in the fits is plotted as a dashed line. Improvement in scatter when fitting for cutout sizes above 20 Petrosian 
Radii is limited, so we use a 20 halflight radii cutout size for all images. Fit parameters seem to have no sensitivity to cutout size in this 
range, suggesting that it may even be possible to use smaller cutouts. 



SDSS pipeline. If we correct the PyMorph estimate of sky 
by the bias noted from the simulations, we find that Py- 
Morph underestimates sky brightness at the 0.1% level while 
the difference between the PyMorph sky estimation and the 
SDSS photometric pipeline is about 0.25%, suggesting that 
the sky values in SDSS are, in fact, overestimated by ap- 
proximately 0.15%. 

Upon further examination of Figure 1151 the bulge pa- 
rameters of the model are more accurately estimated when 
the sky is treated as a free parameter in the fit rather than 
when the sky is fixed at the correct value. The biases dis- 
cussed in Section [3Tl (i. e., underestimate of the Sersic index) 
are partially offset by the biases introduced due to under- 
estimate of the sky. This leads to an apparent improvement 
in the recovery of fit paramters. 



3.5 Testing with real images 

To verify the validity of the simulated background and to 
test the fitting pipeline in clustered environments, we insert 
the simulated galaxies into real SDSS fpC images. The fpC 
images are selected from SDSS DR7 images containing spec- 
troscopic galaxy targets to ensure that any environmental 
changes are minimized. 



In Figures 6(c) 7(c) 8(c) 9(c) and 10(c) we show the 



scatter resulting from fitting the simulated galaxies inserted 
into real fpC images. The scatter and the bias in the fits 
remains unchanged, suggesting that we have properly mod- 
eled the sky background and neighboring sources common 
to an SDSS spectroscopic galaxy. 

Dense environments provide an additional test for our 
pipeline. To select fpC images that contain dense environ- 
ments, we use the GMBCG catalog (Hao et al. 2011). We 
match cluster BCGs with galaxies in our original catalog to 
select fpC images with cluster members including the BCG. 
Our simulated galaxies are then inserted into the image. 
The simulated galaxy in a clustered environment is then 
run through the pipeline. In our previous simulations, intra- 
cluster light and gradients in the sky were not modeled. 
These tests allow us to see what the effects may be. Fig- 
ures [6(d)] [7(d)] [8(d)] [9(d)] and|10(d)|show the results of the 



fits. Once again, the errors remain unchanged, showing that 
no environmental correction is necessary when using the fits 
from the pipeline. 

Placing simulated galaxies near cluster members allows 
us to test for systematic effects in crowded fields. However, 



further examination of BCG galaxies is necessary before 
we are able to properly model them for this purpose. For 
example, the curvature at the bright end observed in the 
size-luminosity relation of early-type galaxies (see Bernardi 
et al. 2012) appears to be due to an increasing incidence 
of BCGs, which define steeper relations than the bulk of 
the early-type population (e.g., Bernardi et al. 2007, 2012). 
However, the curvature could also be due to intracluster light 
(e. g., Bernardi 2009). Further tests for the largest, brightest 
galaxies are needed to explore this possibility. We have not 
presented these tests in this text. 



3.6 Testing HST quality images 

As a final test, we simulated galaxies using settings similar 
to HST data for the F850LP filter band using the Advanced 
Camera for Surveys (ACS). Although similar to testing re- 
covery at higher redshift, the focus is on examining the ef- 
fects of further lowering the S/N. 

Prior to simulation, the sizes and magnitudes of the 
galaxies are rescaled to the apparent size and magnitude of 
the same galaxies at a redshift of z ~ 1.5. However, proper 
k-corrections and size corrections are not applied. The pixel 
size is decreased by a factor of w 10 (from 0.396 to 0.03 
arcsec/pixel). Additional adjustments are made to the zero- 
point and exposure time to properly account for HST-like 
observing conditions. 

For background images, the simulated GEMS tiles are 
used (for further discussion, see Haussler et al. 2007). Fol- 
lowing Haussler et al. (2007), we adopt a hybrid simulation 
method. Standard HST images typically undersample the 
PSF and require the use of drizzling to reconstruct the prop- 
erly sampled image. Nevertheless, we simulate images with 
the final full resolution (0.03 arcsec/pixel). We add Poisson 
noise to the simulated image. Then the simulated galaxy is 
inserted into a simulated image from the GEMS survey. This 
ensures that the galaxies have noise models similar to real 
HST images. 

These simulations allow us to explore lower S/N fits and 
the effects on parameter recovery. The S/N as a function of 
absolute magnitude for the simulations presented in sections 
13. II through 13. 51 and this section are compared in Figure [16] 

The reduction in S/N increases the scatter without bi- 
asing the fits. Examples of the scatter in the total magni- 
tude and halflight radius are plotted in Figure [171 The 68% 
scatter is increased to « 0.05 — 0.10 in apparent magnitude 
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Figure 14. The simulated and recovered apparent magnitude (left), halflight radius (center), and Sersic index (right) for a Ser galaxy fit 
with a Ser model The residuals are plotted as a function of the sky level. Points plotted in black are from fits performed with fixed sky. 
The ovcrplotted points in red are the result of fitting with sky level as a free parameter in the fit. The vertical dashed red lines mark the 
68% scatter of the free sky determination. Our fits are slightly biased low, and this contributes to a small overall bias in fit parameters. 



< 



i.o : 
0.5 ; 
o.o ; 
-0.5 : 
-i.o : 



i i 
i i 



u 
a> 

u 

i— 
rt3 

a 
<i 



4.0 E-l 
2.0 B 
o.o E-"* 

-2.0 Er 
-4.0 EEI 



I I 
I I 



2.0 



< 0.0- 



-2.0 — 






Asky r (%) 



Asky r (%) 



Asky r (%) 



u 

QJ 

u 




Asky r (%) 



Asky r (%) 



u 

QJ 



U 
(13 




Asky r (%) 



c 

< 



2.0 



0.0 




Asky r (%) 



Asky r (%) 




o l 
Asky r (%) 



Figure 15. The simulated and recovered parameters of a simulated SerExp galaxy fit with a SerExp profile. The residuals are plotted as 
a function of the sky level. Points plotted in black are from fits performed with fixed sky. The overplotted points in red are the result of 
fitting with sky as a free parameter of the model. The vertical dashed red lines mark the 68% scatter of the GALFIT sky determination. 
Notice that disk parameters are relatively robust while bulge parameters suffer from incorrect sky estimation. Underestimates of sky 
level have a particularly strong effect on the bulge. 



and « 0.1 arcsec in halflight radii. The dependence of the 
scatter on magnitude parallels the average S/N vs. apparent 
magnitude and reflects the sensitivity to S/N that we have 
referred to throughout the paper. This dependence is also 
apparent in the scatter of the halflight radius when exam- 
ined vs. magnitude. 



4 DISCUSSION 

In the preceding sections we have shown the covariance, bias, 
and scatter in our parameter estimation for the Ser and 
SerExp models. In reality, the effects above will combine 
to yield a total scatter, covariance, and bias that should 
approach those shown in 13.21 Our simulations give us an 



idea of the behavior of the PyMorph pipeline when fitting 
SDSS galaxies as presented in M2013. 

The simulations show that the recovery of global fit- 
ting parameters (total magnitude and halflight radius) in 
the case of SDSS galaxies is remarkably robust, even in the 
case of the SerExp fits. In the case of two-component fits, 
the bulge component is exhibits systematics. However, the 
two component fits still provide valuable insight into the 
structure of galaxies. 

Many galaxies exhibit more complex structure than 
a single or two component structure. Even the case of a 
two component model often oversimplifies galaxy structure. 
Bars, rings, central sources, clumpyness, or asymmetry can- 
not be effectively modeled in our simulations. Because of 
this, we can only determine a lower-bound on the uncer- 
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Figure 17. Simulated and recovered magnitude and halflight radius of a simulated Ser galaxy fit with a Ser profile. Other fits have 
similar behavior. The scatter in recovered magnitude has increased substantially. The shape of the scatter in the magnitude plot mirrors 
the behavior of the S/N. Estimation of the halflight radius exhibits less of an effect, and may even improve, likely due to the increased 
resolution. 
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Figure 16. The distribution of galaxy S/N as a function of ab- 
solute magnitude for the Ser model in red, the DevExp model in 
green, and the SerExp model in blue. The median in each bin 
is shown with symbols. The SDSS galaxies are shown with large 
filled circles, while the redshifted sample is shown with smaller 
filled points. The redshifted sample has S/N lower by a factor of 
w 5. 



tainty in our parameter estimates. However, correcting fits 
using this lower bound represents an improvement in the 
fits of the observed galaxy. Additionally, the galaxies fit in 
M2013 have a median size roughly equivalent to the average 
PSF of SDSS, meaning that, for many galaxies, the resolu- 
tion necessary to resolve this substructure is not present. 

Our tests were performed on r-band data from SDSS. 
The performance of the pipeline can change when observing 
in different bands. This change is primarily dependent on 
the change in the S/N between bands (due to the changing 
brightness of the sky and the color of the galaxy) and on 
the different galactic structures to which neighboring SDSS 
filter bands are sensitive. In M2013, we fit the SDSS g, r, and 
i band data. It is unlikely that the images change drastically 
enough over the wavelength and redshift range observed to 
require additional testing in the g and i bands. However, 
these simulations become an increasingly poor estimate of 
error in bluer bands where the photometry becomes more 
sensitive to star forming regions. These regions tend to be 
clumpier and, therefore, less well represented by a smooth 
profile. Therefore, g band fits may present more scatter than 
the r or i band data. 

The question of whether it is possible to perform a two- 
component fit on data with moderate to low S/N is often 
presented in light of fits similar to these. When a sample 



of galaxy parameters with proper background and S/N are 
fitted, the resulting fits are indeed robust for many parame- 
ters, although specific areas of the parameter space must be 
treated carefully. In particular, the bulge components of a 
two component SerExp fit suffer from systemmatics, partic- 
ularly at larger radii and larger Sersic indexes. However, the 
bulge and disk parameters are recoverable for a wide range 
of simulated model parameters. With this information, we 
can look more seriously at the fits presented in M2013. 

The question of which fit is the "correct" fit for any 
given galaxy is also often asked. To answer this question, 
one can use the reduced \ 2 f° r each fit type to test whether 
a statistically significant change occurs when relaxing the 
parameters of the model (i. e., going from a one-component 
to two-component fit, or allowing the Sersic index to vary) 
or whether improvement in the fit simply results from more 
freedom in the model. This comparison was done by Simard 
et al. (2011). This is certainly a statistically valid question. 
However, the question of which fit is the most appropriate 
fit should not be confused with stating that the "incorrect" 
models do not reveal valid information about the galaxies 
of interest. In Bernardi et al. (2012), we show that the use 
of the all the models jointly can often lift degeneracy and 
reveal important systematic trends in the fitting process. 



5 CONCLUSION 

We presented the simulations used to test fitting of SDSS 
galaxies using PyMorph. Simulations of the Ser and SerExp 
models were presented and examined in many different 
cases. The simulations were generated using the results of 
the fits presented in M2013. We showed that our simulations 
are recoverable in the case of no noise, which demonstrates 
that our simulations are correct. We then showed that we 
can recover the parameters in the case of a simulated back- 
ground and noise representative of the average SDSS image. 

Several individual effects on the fitting were examined. 
We showed that our choice of 20 halflight radii for cutout 
size does not significantly bias our fitting results. In addition, 
we examined the effect of incorrect background estimation, 
which can significantly affect fitting results. 

We also examined the bias created when fitting incor- 
rect models to galaxies. Fitting a two-component Sersic + 
Exponential model to what is really just a single Sersic re- 
sults in a noisier recovery of the input parameters, but these 
are not biased; fitting a single Sersic to what is truly a two- 
component system results in significant biases. These biases, 
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for total magnitude and halflight radius in particular, should 
be useful in correcting other automatic fitting routines. 

These simulations can be analyzed together with the fits 
presented in M20f3 to give a more detailed understanding 
of galaxy structure and formation as presented in Bernardi 
et al. (2012). 
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